令 $f$ 为定义域在实数集上的函数。如果 $\forall x \in \mathbb{R}, f^{\prime\prime}(x) \geq 0$，则称 $f$ 是凸函数。当 $f$ 的变量是向量时，则条件推广为 $f$ 的海森矩阵 $H$ 是正定矩阵（$H \geq 0$）。如果 $\forall x \in \mathbb{R}, f^{\prime\prime}(x) > 0$，则称 $f$ 是严格凸函数（向量函数的情况，对应的条件为 $H$ 是严格正定矩阵，写作 $H>0$）。Jensen不等式的表述如下：
定理：如果 $f$ 是凸函数，令 $X$ 表示某个随机变量，那么 $E[f(X)] \geq f(EX)$
另外，如果 $f$ 是严格凸函数，那么当且仅当以 $1$ 的概率 $X=E[X]$ 时（也即 $X$ 是常数），有 $E[f(X)]=f(EX)$ 。
为了方便记忆，可以想象 $f(x) = x^2$，而 $X$ 是一个服从参数 $0.5$ 的伯努利分布随机变量，取值范围为 $a, b$。那么，$f(EX)=f(\frac{a+b}{2})=\frac{(a+b)^2}{4}$，而 $E[f(x)]=\frac{f(a)+f(b)}{2}=\frac{a^2+b^2}{2}$，显然 $E[f(X)] \geq f(EX)$。
备注：当且仅当 $-f$ 是（严格）凸函数时，$f$ 是（严格）凹函数（即 $f^{\prime\prime}\leq0$ 或 $H \leq 0$）。对于凹函数来说，Jensen不等式的方向改变，$E[f(X)] \leq f(EX)$。
假设对于一个预测问题，有包含 $m$ 个独立样本的训练集 $\{x^{(1)},\cdots,x^{(m)}\}$。我们希望针对数据，估计模型 $p(x,z)$ 的参数，这个模型的对数似然函数由下式给出： $$ \begin{split} \ell(\theta) &= \sum_{i=1}^m log p(x;\theta) \\ &= \sum_{i=1}^m log \sum_z p(x,z;\theta) \end{split} $$ 但是，要显式地找到参数 $\theta$ 的最大似然估计值是十分困难的。这里，$z^{(i)}$ 叫做隐式随机变量。通常如果 $z^{(i)}$ 可以直接观测到，那么最大似然估计法是容易求解的。
在上述情境下，EM算法可以有效地进行最大似然估计。尽管显式地最大化 $\ell(\theta)$ 非常困难，这里的策略是，不断地构造 $\ell$ 的下界（E步骤），然后优化下界（M步骤）。
对每个 $i$，令 $Q_i$ 为 $z$ 的某个分布（$\sum_zQ_i(z)=1，Q_i(z) \geq 0$）。考虑下面的推导： $$ \begin{split} \sum_i log p(x^{(i)};\theta) &= \sum_i log \sum_{z^{(i)}} p(x^{(i)},z^{(i)};\theta) \\ &= \sum_i log \sum_{z^{(i)}} Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \\ &\geq \sum_i \sum_{z^{(i)}} Q_i(z^{(i)}) log \frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \end{split} $$ 最后一步的推导，用到了Jensen不等式。注意到 $f(x)=logx$ 是一个凹函数，因为 $\forall x \in \mathbb{R}^+, f^{\prime\prime}(x)=-\frac{1}{x^2}<0$。而 $$ \sum_{z^{(i)}} Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} = E[\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}] $$ 将期望项中的变量视为 $x$，从而 $$ f(E_{z^{(i)} \sim Q_i}[\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}]) \geq E_{z^{(i)} \sim Q_i}[f(\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})})] $$ 这里的下标 $z^{(i)} \sim Q_i$ 意味着，期望值受限于 $z^{(i)}$ 从分布 $Q_i$ 中取得。
这样，给定任意一组分布 $Q_i$，上面这个不等式都给出了 $\ell(\theta)$ 的下界。但问题是，存在无数种 $Q_i$，我们要怎么选？这样，如果已经有参数 $\theta$ 的某个猜测值，很自然地，我们会根据这个 $\theta$ 至少取到 $\ell(\theta)$ 的下界。（后面我们会发现，这样做使得 $\ell(\theta)$ 在EM的迭代过程中单调上升）
而要使 $\ell(\theta)$ 取得下界，则需要Jensen不等式取等号，而Jensen不等式取等号的充分条件，如第一节所说的，需要随机变量为常数，在当前的情境下，也即 $$ \frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}=c $$ 由于 $Q_i(z^{(i)})$ 是一个概率分布，所以有 $\sum_z Q_i(z^{(i)})=1$，于是，有 $$ \begin{split} Q_i(z^{(i)}) &= \frac{p(x^{(i)},z^{(i)};\theta)}{\sum_z p(x^{(i)},z^{(i)};\theta)} \\ &= \frac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)} \\ &= p(z^{(i)}|x^{(i)};\theta) \end{split} $$ 也就是说，设 $Q_i$ 为设参数为 $\theta$，给定 $x^{(i)}$ 时 $z^{(i)}$ 的后验概率。
现在，在选择了 $Q_i=p(z^{(i)}|x^{(i)};\theta)$ 之后，上面的不等式给了我们想要最大化的 $\ell$ 一个下界，这是E步骤。而在M步骤中，根据上面的不等式，将前面的下界作为一个关于 $\theta$ 的最大值优化问题求解，获得新的参数 $\theta$。重复这两个步骤至收敛，就是EM算法。
怎样确定EM算法会收敛呢？假设 $\theta^{(t)}$ 和 $\theta^{(t+1)}$ 是连续两次EM迭代过程中的参数，我们将证明，$\ell(\theta^{(t)}) \leq \ell(\theta^{(t+1)})$，这也意味着EM算法总是单调提升对数似然。证明的关键，还是在于之前对于 $Q_i$ 的选择。具体来看，迭代起始时，估计的参数是 $\theta^{(t)}$，从而，我们选择 $Q_i^{(t)}(z^{(i)})=p(z^{(i)}|x^{(i)};\theta^{(t)})$，这里会使对应的Jensen不等式取到等号，因此 $$ \ell(\theta^{(t)}) = \sum_i \sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) log\frac{p(z^{(i)}|x^{(i)};\theta^{(t)})}{Q_i^{(t)}(z^{(i)})} $$ 而参数 $\theta^{(t+1)}$ 是通过最大化上面的右式，获得的值，于是有 $$ \begin{split} \ell(\theta^{(t+1)}) &\geq \sum_i \sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) log\frac{p(z^{(i)}|x^{(i)};\theta^{(t)})}{Q_i^{(t+1)}(z^{(i)})} \\ &\geq \sum_i \sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) log\frac{p(z^{(i)}|x^{(i)};\theta^{(t)})}{Q_i^{(t)}(z^{(i)})} \\ &= \ell(\theta^{(t)}) \end{split} $$
备注 如果定义 $$J(Q,\theta)=\sum_i\sum_{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$$ 根据之前的定义，已经知道 $\ell(\theta) \geq J(Q,\theta)$。从而EM算法可以视作是针对 $J$ 的坐标上升法。E步骤根据 $Q$ 求最大值，而M步骤根据 $\theta$ 求最大值。
有了EM算法的理论基础，我们再回过头看高斯混合模型的推导过程。
E步骤较为简单 $$ w_j^{(i)} = Q_i(z^{(i)}=j) = P(z^{(i)}=j|x^{(i)};\phi,\mu,\Sigma) $$ 这里 $Q_i(z^{(i)}=j)$ 表示在 $Q_i$ 的分布前提下，$z^{(i)}$ 取值 $j$ 的概率。
在M步骤，我们需要根据参数 $\phi, \mu, \Sigma$，求下面这个式子的最大值 $$ \begin{split} & \sum_{i=1}^m \sum_{z^{(i)}} Q_i(z^{(i)}) log \frac{p(x^{(i)},z^{(i)};\phi,\mu,\Sigma)}{Q_i(z^{(i)})} \\ &= \sum_{i=1}^m \sum_{j=1}^k Q_i(z^{(i)}=j) log \frac{x^{(i)}|p(z^{(i)}=j;\mu,\Sigma)p(z^{(i)}=j;\phi)}{Q_i(z^{(i)})} \\ &= \sum_{i=1}^m \sum_{j=1}^k w_j^{(i)} log \frac{\frac{1}{(2\pi)^{n/2}|\Sigma_j|^{1/2}}exp(-\frac{1}{2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j)) \cdot \phi_j}{Q_i(z^{(i)})} \end{split} $$
先根据 $\mu_l$ 来计算偏导 $$ \begin{split} &\nabla_{\mu_l}\sum_{i=1}^m \sum_{j=1}^k w_j^{(i)} log \frac{\frac{1}{(2\pi)^{n/2}|\Sigma_j|^{1/2}}exp(-\frac{1}{2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j)) \cdot \phi_j}{Q_i(z^{(i)})} \\ &= -\nabla_{\mu_l}\sum_{i=1}^m \sum_{j=1}^k w_j^{(i)} \frac{1}{2}(x^{(i)}-\mu_j)^T \Sigma_j^{-1}(x^{(i)}-\mu_j) \\ &= \frac{1}{2}\sum_{i=1}^m w_j^{(i)} \nabla_{\mu_l} 2\mu_l^T \Sigma_l^{-1}x^{(i)} - \mu_l^T\Sigma_l^{-1}\mu_l \\ &= \sum_{i=1}^m w_j^{(i)} (\Sigma_l^{-1}x^{(i)} - \Sigma_l^{-1}\mu_l) \end{split} $$ 令上式为零，可以解出 $$ \mu_l = \frac{\sum_{i=1}^m w_j^{(i)} x^{(i)}}{\sum_{i=1}^m w_j^{(i)}} $$
$\phi_j$ 的求解过程，需要利用拉格朗日乘数法。而 $\Sigma_j$ 的求解也是类似的。
In [ ]: